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Abstract 

In this paper numerical calculations are presented for the steady-periodic tempera- 
ture in layered materials and functionally-graded materials to simulate photothermal 
methods for the measurement of thermal properties. No laboratory experiments were 
performed. The temperature is found from a new Green’s function formulation which 
is particularly well-suited to machine calculation. The simulation method is verified 
by comparison with literature data for a layered material. The method is applied to a 
class of two-component functionally-graded materials and results for temperature and 
sensitivity coefficients are presented. An optimality criterion, based on the sensitivity 
coefficients, is used for choosing what experimental conditions will be needed for pho- 
tothermal measurements to determine the spatial distribution of thermal properties. 
This method for optimal experiment design is completely general and may be applied 
to any photothermal technique and to any functionally-graded material. 

KEY WORDS: functionally graded material; Green’s functions; optimal experiment; 
photothermal; thermal properties 


1 Introduction 

Functionally-graded (FG) materials are being studied as possible components of aero- 
space thermal protection systems. These materials include composites with epoxy and 
metal matrices, metal foams, or any structure with properties designed to vary with 
position. In the future when FG materials are specified as part of a vehicle program, 
part of the procurement process will involve certification that the material meets the 
specifications. 

To date there has been little research on accurate thermal characterization of FG 
materials. The present research is intended to close this gap in the procurement cycle 
by investigating photothermal methods for non-destructive and accurate measurement 
of thermal properties in FG materials. In this paper only numerical simulations are 
presented and no laboratory experiments were performed. 
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2 Mechanical Engineering Dept., University of Nebraska-Lincoln, Lincoln, NE 68588-0656 
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A review of the pertinent literature is given next in three areas: computer simu- 
lation of FG materials; heat transfer theory for photothermal methods; and, optimal 
experiment design. 

Several researchers have found analytical solutions for the thermal response by 
representing a FG material as composed of multiple layers each with different, spa- 
tially uniform, thermal properties [1,2,3]. Other studies investigated an exponential- 
function variation of thermal properties along one spatial direction [4,5]. Another used 
Galerkin’s method to find temperature in materials with arbitrary property distribu- 
tions [6]. The primary motivation for these studies has been to determine temperature 
for the purpose of finding thermal stresses, or to find the distribution of thermal prop- 
erties that optimizes the thermal stresses [7,8]. 

One research group has reported transient-heating experiments to measure thermal 
properties in a FG material [9,10]. This group studied a FG material containing an 
exponentially-varying spatial distribution of thermal properties. Their data analysis 
combines a single temperature datum with their transient theory to provide a sin- 
gle value for the parameter describing the spatial distribution of thermal properties. 
Although simple in concept, this approach is sensitive to measurement noise. 

In the area of photothermal measurements, there are several pertinent publica- 
tions. A diverse collection of thermal-wave Green’s functions and temperature solu- 
tions has been published recently in book form [11]. Primarily homogeneous materials 
are treated, and layered materials are included by defining a global Green’s function 
that embodies the effects of several layers in the material. Because the complexity of 
the layered-body Green’s function increases rapidly as layers are added, no more than 
3 layers are discussed. 

Theory for many-layered bodies to laser heating has been studied previously by the 
author [12,13]. The volumetric heating is treated exactly from the optical absorption 
properties of all layers. The multi-layer body is treated efficiently by the use of local 
Green’s functions which are found first in the time domain and are then transformed 
into the frequency domain. Each layer is linked to adjacent layers with appropriate 
interface conditions. 

Theory for the photoacoustic response of a layered solid has been recently studied 
for which the temperature in each layer is linked with adjacent layers by interface 
conditions [14]. The optical absorption in each solid layer is described by an exponen- 
tial distribution and an absorption coefficient. The photoacoustic response is found 
from both thermal effects and mechanical effects in the gas, however thermal effects 
predominate for solid materials. The method is used for analysis of experimental data 
in materials with 2 and 3 layers. 

In the area of optimal experiment design, parameter estimation has been used for 
obtaining thermal properties from transient experiments for many years [15]. In these 
methods the desired parameters are found by non-linear regression between the experi- 
mental data (temperatures in this case) and a computational model of the experiment. 
Parameter est im ation concepts have recently been applied to optimal experiment de- 
sign for thermal characterization of uniform materials [16] and for materials with 
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temperature- varying properties [17]. The author has previously studied optimal ex- 
periment design for low-conductivity FG materials [18]. The simulated experiments 
involved time-series data collected from one or more temperature sensors, and the 
data analysis is carried out in the time domain. An optimality criterion was used to 
find the best experimental conditions for simultaneous estimation of several thermal 
properties. Results of simulations show that for FG materials with spatially-varying 
conductivity, it is better to heat the sample from the low-conductivity side. 

In the present paper, FG materials are simulated with a large number of intercon- 
nected layers. The heat transfer theory draws upon the author’s previous work with 
Green’s functions, but here the Green’s functions are given directly in the frequency 
domain in the form of algebraic expressions, not infinite-series expressions, that are 
numerically well behaved under all conditions. Likewise the temperature expressions 
found from these Green’s functions are numerically well behaved. The Green’s func- 
tions are given for a variety of boundary conditions; previously only specified-flux 
boundaries were treated. To the author’s knowledge this paper describes the first 
application of optimal experiment design methods to frequency-domain analysis ap- 
propriate for photothermal experiments. This paper is divided into several sections, 
as follows: the temperature in one layer; the Green’s functions; the temperature in a 
multi-layer material; the design of optimal experiments; results for a layered material; 
results for a FG material; and, a brief summary. 


2 Temperature in one layer 

Since the photothermal applications of interest involve periodic heating by a laser, the 
solution is sought in Fourier-transform space, and the solution is interpreted as the 
steady-periodic response at a single frequency ui. For a discussion of this point see 
[19]. Consider the heat conduction equation in Fourier transform space in one layer: 

§--f r = 0<X<L (1) 

dT 

ki~ VhiT — fi{u)-, at boundaries i = 1,2 (2) 

orii 

Here T is Fourier-space temperature (K s), a is thermal diffusivity (m 2 /s), k is thermal 
conductivity (W/m/K), g is volume heating (W s/m 3 ) deposited by a laser, and /* is a 
specified boundary condition. Index i = 1,2 represents the boundaries at the limiting 
values of coordinate x. The boundary condition may be one of three types at each 
boundary: for type 1 fo is a specified temperature (k{ = 0 and hi = 1); for type 2 /, is 
a specified heat flux (hi — 0); and, type 3 represents a convection condition where h t 
is a constant- with-time heat transfer coefficient (or contact conductance). 

The temperature will be found with the Fourier-space Green’s function, defined by 
the following equations: 

~ ° 2 G = -~&(x-x') (3) 

ox z a 
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(4) 


ki~ 1- hiG = 0; i — 1,2 

OTli 

Here a 2 = ju/a and 6(x-x') is the Dirac delta function. The coefficient 1/a preceding 
the delta function in Eq. (3) provides the frequency-domain Green’s function with 
units of seconds/meters. This is consistent with our earlier work with time-domain 
Green’s functions. 

Assume for the moment that the Green’s function G is known, then the steady- 
periodic temperature is given by the following integral equation (see [20], p. 40-43): 


T{x,uj) 


J g(x', cj)G(x, x',u)dx' (for volume heating) 


+afi(u>) x 


dG/dn'(x,Xi,oj) 

\G{x,Xi,u)) 


(type 1 only) 
(type 2 or 3) 


i = 1,2 


( 5 ) 


3 Green’s Function 


The Green’s function (GF) that satisfies Eqs. (3) and (4) is given by 


G(x,x',uj) = 


+ 


S 2 (ST e~ <r ( 2L - |x - x 'l ) + 5 1 + e~ <T(2L_x_a: ' ) ) 
2aa{S/St - SfS^eF^ 1 ) 
Sj(5fe-^ | *"* ,|) + Si e -*(*+*'>) 

2aa{S± S} - 5f5 2 -e-^) 


( 6 ) 


where the subscripts 1 and 2 represent the two boundaries at the smallest and largest 
x-values, respectively. Coefficients S m and depend on the boundary conditions on 
side M and are given by 


‘S’m 




1 

kc T + hf4 


if side M is type 0, type 1, or type 2 
if side M is type 3 


0 if side M is type 0 

—1 if side M is type 1 

1 if side M is type 2 

k ka — h M if side M is type 3 


A boundary of type 0 designates a far-away boundary, as in a semi-infinite body. The 
derivation of the Fourier-space GF in Eq. (6) parallels that for steady-state GF given 
elsewhere [21]; however in the present work a is complex. 

This form of the GF is particularly well-behaved for machine computation, and 
most importantly, the temperature expressions based on these GF are similarly well- 
behaved for any layer thickness and for any frequency. This is a key contribution of this 
paper, in sharp contrast with previously reported difficulties in evaluating numerical 
values from exact solutions. For example, numerical overflow can occur with other 
formulations in thermally thick layers [14]. In a time-domain study, only short-time 
results were included due to numerical difficulties associated with longer times [3]. 
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The GF expression given in Eq. (6) covers a number of boundary condition com- 
binations, and a numbering system is used to distinguish among them. Designation 
XIJ is used to identify the GF for heat transfer in a layer with boundary condition 
of type I = 1, 2, or 3 at x = 0 and with boundary condition of type J = 1, 2, or 3 
at x = L. For example, designation X12 represents the GF with type 1 boundary at 
x = 0 and type 2 boundary at x — L. Designation XIO is used to identify the GF for 
a semi-infinite region with a boundary of type I = 1, 2, or 3 at a; = 0. 

4 Temperature in Layered Materials 

In this section the temperature caused by absorption of laser energy will be found 
in a domain consisting of non-absorbing air, N solid layers, and a substrate. At 
the interfaces, let q nm represent the heat flux leaving layer n and entering layer m. 
Applying Eq. (5), the interface temperature in the air is: 

T o (0, w) = — ^Go(0, 0, u)q w (7) 

Kq 

In layer v, i — 1,2, N : the interface temperatures axe: 

Ti(0,u>) = -j^Gi(0,0,u)qi-i t i + -j^Gi(0,Li,u})q i+ i )i + Bi(0) (8) 

r V| rCj 

Ti(Li,u ;) = -£-Gi(Li,0,u)qi-ij + ^-Gi(Li,Li,u)qi + i t i + Bi(Li) ( 9 ) 

/Cj Ki 

In the substrate the temperature at the interface is: 

Tv +1 (0,u;)= Gn+i (0, 0, ^)gw,N + i + -Efyr+i (0) (10) 

«w+i 

In the above expression, symbol Bi has been used for the volume-heating integral term 
from Eq. (5), specifically, 

Bi{x) = -r [ g{x',u) Gi(x,x',u}) dx' (11) 

Ki Jx* 

Here g is the laser energy absorbed in the layer per unit volume; this can be determined 
without approximation from the optical properties of the layers [12], 

In the above temperature expressions, all of the interface heat fluxes are initially 
unknown. The heat flux leaving one layer enters the adjacent layer, = —q^i-x 

and the temperature difference between adjacent layers is related to heat flux through 
a contact resistance at each interface: 

Qi-i,i Ri = Ti(0, cj) — Ti-i(Li-i, u); i — 1, 2, . . . , N + 1 (12) 

The contact resistance Ri describes the size of the temperature jump across the inter- 
face. Next Eqs. (7-10) are combined with Eq. (12) to eliminate temperature. The 


result is a set of N + 1 linear algebraic equations for the unknown heat fluxes, which 
may be stated in matrix form: 

Wo U i R\ — Vi 0 ... 0 

—Vi Ui + U 2 + R 2 — V 2 ••• 0 

0 —V 2 U 2 + U 3 4- ••• 0 

-V N 

0 0 ... —Vff Un + W/v+i + Rn+1 



010 


' ^(0) 


921 


B 2 ( 0 ) - Bi(L\) 

X 

932 

— 

B s (0) — B2(L 2 ) 


. Qn+i,n . 


. Bn + i(0) - Bh(Ln) 


Symbols Wi, Uj, and V* used in the above expression are given below: 

w, = 2t G( (0,0,w) 

U, = %G l (0,a,u)=%G i (L l ,L„u) 

V = %-G i (0,L i ,u>) = ^G i (L i ,0,u) (14) 

For any multilayered system, it is now possible to calculate the heat fluxes (q t j) through 
all interfaces in the system. The above result is exact, and Cramer’s rule may be used 
to solve for the q's for a sample composed of one or two layers. For a sample with two 
or more layers, a numerical matrix solution is best. Once the heat fluxes are found, 
the temperature at any interface is given by Eq. (8 - 10), or the temperature within 
any layer may be found with Eq. (5). 

Several different GF may be used in the above matrix equation. The non-absorbing 
gas (region 0) is a semi-infinite region so the GF needed is number X20. For layers 
i = 1, 2, . . ., N the GF needed are type X22 (specified heat flux). The GF for 
the substrate depends on the heat transfer environment there. For example, a thick 
substrate could be described by GF number X20, or, a substrate in imperfect contact 
with a cold plate could be described by GF number X23. 

5 Optimal Experiment Design 

Sensitivity coefficients are central to the design of optimal experiments for thermal 
property evaluation. The sensitivity coefficients are defined by 

Xrt 0 = b (15) 

which is the sensitivity for the kih parameter, the jth temperature sensor, and the «th 
frequency. Parameters 6* may include conductivity, specific heat, density, etc. In this 
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research the sensitivity matrix has been computed from the real-valued amplitude and 
real- valued phase of the complex temperature at r frequencies, for which the distinct 
amplitude and phase values are treated as 2 r measurements from each sensor. 

The sensitivity coefficients were computed with a finite-difference procedure to 
approximate the derivative, as follows: 


eb k 


(16) 


Here Ty is the temperature at the ith frequency for the jth sensor. The value of 
e = 0.001 was found to give well-behaved values for X. 

There are two specific requirements that the sensitivity coefficients must satisfy. 
First, the sensitivity coefficients should be as large as possible. Second, when two or 
more parameters are to be measured in the same experiment, the sensitivity coefficients 
must be linearly independent. A formal procedure to quantify these two requirements 
may be constructed if the sensitivity coefficients are assembled into a sensitivity matrix 
X, which is then multiplied by its transpose, given formally by X T X, of size [p x p]. 
The optimality criterion is the (normalized) determinant of matrix X T X, given by 
[16] 

D = 77(t^ det(xTx> (n) 

Note that optimality criterion D is normalized by the maximum temperature rise 
(squared), the number of sensors s, and the number of frequencies r. This is important 
so that D may be used to compare different experiments. When D is large, the 
sensitivity coefficients will be large and linearly independent [17]. 


6 Results for a Layered Material 

In this section the techniques of simulation and optimal experiment design are applied 
to a SiC >2 film on a Si substrate. This material has been studied elsewhere by a pho- 
toacoustic technique over a frequency range from 2 to 20 kHz [14]. An opaque coating 
of Ni of 20 nm thickness is added to the sample to improve the optical absorption. 
Published results for the thermal conductivity of the SiC >2 film and the contact resis- 
tance between the film and the Si substrate are: k = 1.52 W/m/K and R < 10 -9 K 
m 2 /W. 

Using the present methods, the computed phase of the temperature agreed with 
published values, thus verifying the present approach. The sensitivity coefficients for 
the phase of the surface temperature to variations in the conductivity of the the SiC >2 
layer, k , and the contact resistance with the Si substrate, R, are plotted in Fig. 1. 
The property values used are given in Table 1. The sensitivity to conductivity k has 
a maximum (negative) value at about 6000 Hz. The shape of the R-sensitivity curves 
and the location of the largest value depends strongly on the R value. For R = 10 -8 
and R = 10“ 7 , the sensitivity curves are relatively small, with a broad peak near 3 
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Figure 1: Phase sensitivity to thermal conductivity and to layer-substrate contact 
resistance for a layer of SiC >2 on silicon at various frequencies. 

kHz. For R = 10 -6 however, the sensitivity values are large and positive at small 
frequencies and the curve slopes down to negative values as frequency increases. Some 
conclusions for R < 10 -7 are that the sensitivity coefficient is small, but the frequency 
range considered (0.5 - 20 kHz) captures the peak sensitivity. For R = 10 -6 the largest 
sensitivity may lie outside this frequency range. 

Table 1. Properties used in calculations of a three-layer solid which is 
heated at the Ni surface and exposed to a layer of air on either side. 


layer 

i 

d (fi m) 

k (W/m/K) 

a (m 2 /s) 

R (K m 2 /W) 

air 

0 

- 

0.0263 

2.25E-05 

- 

Ni 

1 

0.07 

80 

1.98E-05 

0 . 

Si0 2 

2 

0.485 

1.52 

9.09E-07 

0 . 

Si 

3 

382 

151 

9.09E-04 

varies 

air 

4 

- 

0.0263 

2.25E-05 

- 


Next the sensitivities for k and R will be examined together. Optimality criterion 
D provides a numerical measure of the extent to which the sensitivity coefficients are 
both large and linearly independent, and at which frequencies this occurs. Figure 
2 shows the optimality criterion for the layered material for both conductivity and 
contact resistance considered together. The highest curve is for k and R = 10 -6 , and 
the peak for this curve is about 4 kHz. For the R = 10 -7 and R = 10 -8 curves the 
peak occurs around 8 kHz. Figure 2 indicates that the best experiment to measure 
both k and R when R = 10~ 6 includes data at 4 kHz. For smaller i?-values the most 
important frequency is 8 kHz, however because D is smaller the analysis of the data 
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Figure 2: Optimality criterion D for simultaneous estimation of thermal conductivity 
and contact resistance, versus frequency, for a layer of Si 02 on silicon. 

may be more difficult. It is important to note that each point on Fig. 2 represents 
a value for D computed from a range of data extending from 20 kHz down to that 
point; data is added from high frequency to low frequency. To repeat, the D-values 
indicate what frequency range is needed for optimal estimation of both k and R from 
experimental data. 


7 Results for a Functionally Graded Material 

In this section the methods for experimental design are applied to a two-phase ce- 
ramic/ceramic material with graded volume fraction of the components. The volume 
fraction profile is assumed to have the form: 

VU*) = 1 - ( x/LY ; V 2 (x) = ( x/Lf (18) 

where Vi and V 2 are the volume fraction of the components. At location x = 0 the 
material is pure component 1 and at x = L the material is pure component 2. The 
particular material considered is composed of TiC and SiC and the thermal properties 
are given in Table 2. 

Table 2. Material properties of TiC and SiC. 



material 

k (W/m/K) 

9 (kg/m 3 ) 

c (J/kg/K) 

1 

TiC 

20 

4900 

700 

2 

SiC 

60 

3200 

1000 
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Figure 3: Thermal conductivity distribution in a functionally-graded material for sev- 
eral values of distribution parameter p. 


The thermal conductivity of the material is given by 


k(x) = ki 


. V2(x)(k 2 -k l ) ' 

ki + ih-k^V l (x)/3 


(19) 


where subscripts 1 and 2 stand for the properties of components 1 and 2, respectively. 
Figure 3 shows the spatial distribution of conductivity for the TiC/SiC material for 
several values of exponent p. The mass density and specific heat are determined by 
the rule of mixtures: 


p{x) = V\{x)p\ -I- V 2 (x)p 2 (20) 

c{x) = Vi(*)c! + V 2 {x)c 2 (21) 

The thermal-stress behavior of this material has previously been studied [3]. 

In photothermal methods, the sample is heated by a periodically modulated laser 
beam, and the surface temperature (or a subsequent acoustic signal) is measured 
at the modulation frequency. The temperatures reported below are computed with 
the layered GF method with 50 layers used to simulate the spatial variation in the 
sample. Both surfaces of the sample are exposed to air. The surface temperature is 
shown in Fig. 4a (amplitude) and Fig. 4b (phase) versus dimensionless frequency. 
The frequency is normalized as /* = fL 2 /a av where L is the material thickness and 
ot a v = (<*1 + oc 2 )/ 2 , the average of the component values. In Fig. 4 the sample is 
heated at x = 0 and the temperature is reported at x = 0, the low-fc side. This 
heating condition provides slightly higher temperature response than for heating on 
the high-A; side (at x — L). Note that in Fig. 4a the amplitude curves are monotonic 
with few distinguishing features. In contrast the phase has a distinctive maximum 
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Figure 4: Surface temperature (a) amplitude, and (b) phase, versus dimensionless 
frequency for a functionally-graded material heated periodically at the surface. 


for each p value, which supports the experimental observation that phase is more 
important than amplitude for photothermal measurement of thermal properties. 

Figure 5a shows the sensitivity of the phase of temperature to k 2 , the thermal 
conductivity of component 2, for three values of spatial distribution parameter p. Note 
that the sensitivity is largest for p = 1, the linear ^-distribution, and is small for other 
distributions. In Fig. 5a the peak sensitivities in a range of dimensionless frequencies 
below /* = 1. This range of frequencies represents thermal waves (generated by 
the periodic heating) that penetrate all the way through the sample thickness. The 
sensitivity to k\, the component 1 thermal conductivity, is not shown because it is 
similar in size and shape to the k 2 sensitivities (but values are negative). 

Figure 5b shows the sensitivity of the phase of the temperature to spatial distri- 
bution parameter p. The curve for p = 1 has a large positive peak at about /* = 0.2 
and a negative peak at about /* = 1.5. The largest sensitivities are again located in 
the range /* < 1. 

The sensitivity to a single parameter is useful for determining a single parameter 
from an experiment. When two or more parameters are to be determined, optimality 
criterion D is instructive. Figure 6a shows optimality criterion D for both conductivity 
k 2 and exponent p. Both amplitude and phase information was used to compute these 
values. The largest curve is for p = 1 which indicates that this spatial distribution 
will provide better estimates of thermal properties k 2 and p from an experiment than 
for other p-values. The location of the peak for each curve indicates which frequencies 
should be included in an experiment. Since the D-values shown were computed from 
high-to-low frequencies, the peak indicates the lowest frequency of data that is needed 
for optimal data analysis. 






Figure 5: Sensitivity of the phase of temperature to (a) conductivity k 2 , and (b) 
exponent p, for the functionally-graded material heated periodically at the surface. 


Optimality criterion D was also investigated for other combinations of parameters. 
The D-v alues computed for p and k\ together, not shown, were identical to Fig. 6a, 
which is expected from the similarity in the shape of the sensitivities to k\ and k 2 
mentioned earlier. The O- values computed for k\ and k 2 together are shown in Fig. 
6b. The important feature of this figure is the vanishingly small values, five orders of 
magnitude smaller than Fig. 6a values, which indicates that k\ and k 2 should not be 
sought simultaneously. As mentioned before the sensitivity coefficients for ki and k 2 
have a similar shape, so the fact that their combined D - value is near zero reinforces the 
idea that the sensitivities for k\ and k 2 are not linearly independent. Finally, D-values 
were also computed for p, k\, and k 2 considered simultaneously (not shown), and these 
values were predictably near zero because of the dependence of ki and k 2 . 

8 Summary 

This paper investigates photothermal methods for thermal characterization of func- 
tionally graded materials through numerical modeling and experiment design. No 
laboratory experiments are reported. A new formulation is given for the temperature 
response of layered materials to periodic heating, based on the method of Green’s 
functions, which is numerically better behaved than previous work. The method has 
been applied to a SiC >2 layer on silicon and compared to literature values to validate 
the method. Optimality criterion D indicates which frequency range of experimental 
data should provide the best possible estimates of the layer conductivity and contact 
conductance. 

The new methods have also been applied to a two-component functionally-graded 
material with a power-law distribution of thermal properties. The largest temperature 
response is found by heating the sample on the low-k side, and the phase of the tem- 
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Figure 6: Optimality criterion D for simultaneous estimation of (a) conductivity k 2 
distribution parameter p, and (b) conductivities k\ and k 2 , for the functionally-graded 
material. 


perature is particularly important for estimation of thermal properties. Component 
conductivities k\ and k 2 have similar-shaped sensitivity coefficients, and consequently 
both cannot be estimated simultaneously from experimental data. The most impor- 
tant parameter is p which describes the spatial distribution of thermal properties in 
the functionally-graded material. Values for optimality criterion D indicate that val- 
ues for p may be found simultaneously with one of the conductivities, but not both. 
Dimensionless frequencies in the range /* < 1 are important for measurement of 
spatial-distribution parameter p. The magnitude of the optimality criterion D also 
suggests that it will be easier to estimate parameters for p rj 1 (near-linear spatial 
variation) compared to other values of p. 
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